---
title: "ATEPS"
author: "Edmund Kelly"
date: "2023-05-23"
output: html_document
---

```{r setup}

setwd("C:/Users/edmun/OneDrive - Nexus365/ATEPS Data/")

```

```{r Packages, echo = FALSE} 
library(OpenMx)
library(umx)
library(haven)
library(tidyverse)
library(stargazer)
library(sjPlot)
library(ggplot2)
library(broom)
library(RColorBrewer)
library(psych)
library(modelsummary)
library(fixest)
library(psychTools)
set.seed(1234)

```

```{r Loading data}

ATEPS_orig = read_dta("C:/Users/edmun/OneDrive - Nexus365/ATEPS Data/ADA ATEPS Version 6 Feb 2023.dta")

```

```{r Pre-Processing}
 
#Filtering

ATEPS = ATEPS_orig %>% 
  dplyr::select(tra_zygosity, gender, age, male, person_id, twin_id, num_twins, MZ, DZ, twin_years, trust_sent, starts_with("trust_return"),  
         starts_with("pol"), trust_avg_return, peas_pod_q1, peas_pod_q2, first_res, complete_pair, trust_stated)

#Removing missing values and incomplete pairs
ATEPS[ATEPS == 999] <- NA
ATEPS[ATEPS == 99999999] <- NA


#Political trust

ATEPS$pol_job_well = ATEPS$pol_job_well %>% 
  case_match(3 ~ 1, 
             2 ~ 0.5, 
             1 ~ 0)

ATEPS$pol_dont_care = ATEPS$pol_dont_care %>% 
  case_match(3 ~ 0, 
         2 ~ 0.5, 
         1 ~ 1)
ATEPS$pol_misuse = ATEPS$pol_misuse %>% 
  case_match(3 ~ 0, 
         2 ~ 0.5, 
         1 ~ 1)

ATEPS$pol_power = ATEPS$pol_power %>% 
  case_match(3 ~ 0, 
         2 ~ 0.5, 
         1 ~ 1)

ATEPS$trustpol <- as.numeric(ATEPS$pol_dont_care + ATEPS$pol_job_well + ATEPS$pol_misuse + ATEPS$pol_power)/4

ATEPS$trustpol_scaled = as.numeric(scale(ATEPS$trustpol))

ATEPS <- subset(ATEPS, complete_pair == 1) 

#Behavioral trust

ATEPS$trust_sent = 1 - (ATEPS$trust_sent)/11

ATEPS$trust_sent_scaled = as.numeric(scale(ATEPS$trust_sent))

ATEPS$trust_stated_scaled = as.numeric(scale(ATEPS$trust_stated))

#Equal environment variables

ATEPS$peas_pod_q1 = case_when(
  ATEPS$peas_pod_q1 == 1 ~ 3, 
  ATEPS$peas_pod_q1 == 2 ~ 1, 
  ATEPS$peas_pod_q1 == 3 ~ 2
)

#Creating twin ID variable
ATEPS$FAMID = ATEPS$twin_id
ATEPS$TWINID = as.factor(ATEPS$first_res)
ATEPS <- ATEPS %>% mutate(TWINID = 
                          dplyr::recode(TWINID, 
                                 "0" = 1,
                                 "1" = 2))

ATEPS$ZYG = as.factor(case_when(
  ATEPS$MZ == "1" ~ 1, 
  ATEPS$DZ == "1" ~ 2
))

#Stacking data by twin pairs
twin1 = ATEPS %>% 
  filter(TWINID == 1)

twin2 = ATEPS %>% 
  filter(TWINID == 2)

ATEPS_final = merge(twin1, twin2, by.x = "FAMID", by.y = "FAMID", suffixes = c("_T1", "_T2")) 

#Creating MZ and DZ datasets
ATEPS_MZ <- subset(ATEPS_final, ZYG_T1 == 1)
ATEPS_DZ <- subset(ATEPS_final, ZYG_T1 == 2)

ATEPS_orig_MZ = ATEPS_orig %>% 
  filter(MZ == 1 & complete_pair == 1) %>% 
  mutate(pol_job_well = as.numeric(pol_job_well), 
  pol_dont_care = 3 - pol_dont_care, 
  pol_misuse = 3 - pol_misuse, 
  pol_power = 3 - pol_power, 
  trustpol = pol_job_well + pol_dont_care + pol_misuse + pol_power, 
  trustpol_scaled = scale(trustpol), 
  pol_job_well_scaled = scale(pol_job_well), 
  pol_misuse_scaled = scale(pol_misuse), 
  trust_sent_scaled = scale(trust_sent), 
  trust_stated_scaled = scale(trust_stated),
  pol_power_scaled = scale(pol_power), 
  pol_dont_care_scaled = scale(pol_dont_care)
  )

ATEPS_orig_DZ = ATEPS_orig %>% 
  filter(MZ == 0 & complete_pair == 1) %>% 
  mutate(non_mainstream = case_when(
    vote_party == 1 ~ 0, 
    vote_party == 4 ~ 0, 
    TRUE ~ 1
  ), 
  pol_job_well = as.numeric(pol_job_well), 
  pol_dont_care = 3 - pol_dont_care, 
  pol_misuse = 3 - pol_misuse, 
  pol_power = 3 - pol_power, 
  trustpol = pol_job_well + pol_dont_care + pol_misuse + pol_power, 
  trustpol_scaled = scale(trustpol), 
  pol_job_well_scaled = scale(pol_job_well), 
  pol_misuse_scaled = scale(pol_misuse), 
  pol_power_scaled = scale(pol_power), 
  pol_dont_care_scaled = scale(pol_dont_care), 
  trust_sent_scaled = scale(trust_sent), 
  trust_stated_scaled = scale(trust_stated)
  )

```

```{r Appendix 2: Descriptive statistics}

ATEPS_summary_MZ = ATEPS %>% 
  filter(tra_zygosity == 'MZ') %>% 
  select(age, male, starts_with("pol"), starts_with("peas"), complete_pair, trust_stated, trustpol, trust_sent, trust_avg_return) %>% 
  as.data.frame() 

stargazer(ATEPS_summary_MZ, iqr = T, median = T, out = "ATEPS_heritability_summary_MZ.htm")

ATEPS %>% 
  filter(tra_zygosity == 'DZ') %>% 
  select(age, male, starts_with("pol"), starts_with("peas"), complete_pair, trust_stated, trustpol, trust_sent, trust_avg_return, first_res) %>%
  as.data.frame() %>% 
  stargazer(iqr = T, median = T, out = "ATEPS_heritability_summary_DZ.html")

```

```{r Appendix 1: Factor analysis trust}

ATEPS_FA = ATEPS_orig %>% 
  mutate(pol_job_well = 
           3 - pol_job_well)

ATEPS_alph = ATEPS_FA %>% 
  select("pol_power", "pol_misuse", "pol_dont_care", "pol_job_well") %>% 
  psych::alpha() 

ATEPS_fact = ATEPS_FA %>% 
  select("pol_power", "pol_misuse", "pol_dont_care", "pol_job_well") %>% 
  na.omit() %>% 
  fa(rotate = "varimax")

fa2latex(ATEPS_fact)

ATEPS_alph

```

```{r Appendix 2: Trust in politicians ACE} 

set.seed(1234)
ACE <- umxACE("ACE", selDVs = "trustpol_scaled", selCovs = c("male", "age"), sep = "_T", mzData = ATEPS_MZ, dzData = ATEPS_DZ, tryHard = "ordinal", optimizer = "CSOLNP")
ACE <- umxConfint(ACE, run = TRUE, level = .95, parm = c("top.A_std[1,1]", "top.C_std[1,1]", "top.E_std[1,1]"),
                  
                  wipeExistingRequests = F, optimizer = "CSOLNP")
summary(ACE, verbose = TRUE)

ACEc <- umxACE("ACE", selDVs = "trustpol_scaled", selCovs = c("male", "age", "peas_pod_q1", "peas_pod_q2"), sep = "_T", mzData = ATEPS_MZ, dzData = ATEPS_DZ, tryHard = "ordinal", optimizer = "CSOLNP")
set.seed(1234)
ACEc <- umxConfint(ACEc, run = TRUE, level = .95, parm = c("top.A_std[1,1]", "top.C_std[1,1]", "top.E_std[1,1]"), 
                  wipeExistingRequests = F, optimizer = "CSOLNP")
summary(ACEc, verbose = TRUE)

AE <- umxACE("AE", selDVs = "trustpol_scaled", selCovs = c("male", "age"), sep = "_T", mzData = ATEPS_MZ, dzData = ATEPS_DZ, optimizer = "CSOLNP")
set.seed(1234)
AE = umxModify(AE, update = c("c_r1c1"), comparison = TRUE)
AE <- umxConfint(AE, run = TRUE, level = .95, parm = c("top.A_std[1,1]", "top.E_std[1,1]"), optimizer = "CSOLNP", 
                  wipeExistingRequests = F)
summary(AE, verbose = TRUE)

CE <- umxACE("CE", selDVs = "trustpol_scaled", selCovs = c("male", "age"), sep = "_T", mzData = ATEPS_MZ, dzData = ATEPS_DZ, optimizer = "CSOLNP")
set.seed(1234)
CE = umxModify(CE, update = c("a_r1c1"), comparison = TRUE)
CE <- umxConfint(CE, run = TRUE, level = .95, parm = c("top.C_std[1,1]", "top.E_std[1,1]"), optimizer = "CSOLNP", 
                  wipeExistingRequests = F)
summary(CE, verbose = TRUE)

```

```{r Appendix 2: Interpersonal trust ACE} 
ACE <- umxACE("ACE", selDVs = "trust_stated_scaled", selCovs = c("male", "age"), sep = "_T", mzData = ATEPS_MZ, dzData = ATEPS_DZ, tryHard = "ordinal", optimizer = "SLSQP")
ACE <- umxConfint(ACE, run = TRUE, level = .95, parm = c("top.A_std[1,1]", "top.C_std[1,1]", "top.E_std[1,1]"),
                  wipeExistingRequests = F,  optimizer = "SLSQP")
summary(ACE, verbose = TRUE)

ACEc <- umxACE("ACE", selDVs = "trust_stated_scaled", selCovs = c("male", "age", "peas_pod_q1", "peas_pod_q2"), sep = "_T", mzData = ATEPS_MZ, dzData = ATEPS_DZ, optimizer = "SLSQP")
set.seed(1234)
ACEc <- umxConfint(ACEc, run = TRUE, level = .95, parm = c("top.A_std[1,1]", "top.C_std[1,1]", "top.E_std[1,1]"),
                  wipeExistingRequests = F,  optimizer = "SLSQP")
summary(ACEc, verbose = TRUE)

AE <- umxACE("AE", selDVs = "trust_stated_scaled", selCovs = c("male", "age"), sep = "_T", mzData = ATEPS_MZ, dzData = ATEPS_DZ, optimizer = "CSOLNP")
set.seed(1234)
AE = umxModify(AE, update = c("c_r1c1"), comparison = TRUE)
AE <- umxConfint(AE, run = TRUE, level = .95, parm = c("top.A_std[1,1]", "top.E_std[1,1]"), 
                  wipeExistingRequests = F, optimizer = "CSOLNP")
summary(AE, verbose = TRUE)

CE <- umxACE("CE", selDVs = "trust_stated_scaled", selCovs = c("male", "age"), sep = "_T", mzData = ATEPS_MZ, dzData = ATEPS_DZ, optimizer = "SLSQP")
set.seed(1234)
CE = umxModify(CE, update = c("a_r1c1"), comparison = TRUE)
CE <- umxConfint(CE, run = TRUE, level = .95, parm = c("top.C_std[1,1]", "top.E_std[1,1]"), 
                  wipeExistingRequests = F, optimizer = "SLSQP")
summary(CE, verbose = TRUE)

```

```{r Appendix 2: Behavioral trust ACE}

ACE <- umxACE("ACE", selDVs = "trust_sent_scaled", 
              selCovs = c("male", "age"), 
              sep = "_T", mzData = ATEPS_MZ, dzData = ATEPS_DZ, tryHard = "yes", optimizer = "CSOLNP")
ACE <- umxConfint(ACE, run = TRUE, level = .95, parm = c("top.A_std[1,1]", "top.C_std[1,1]", "top.E_std[1,1]"),
                  wipeExistingRequests = F,  optimizer = "SLSQP")
summary(ACE, verbose = TRUE)

ACEc <- umxACE("ACE", selDVs = "trust_sent_scaled", 
               selCovs = c("male", "age", "peas_pod_q1", "peas_pod_q2"), 
               sep = "_T", mzData = ATEPS_MZ, tryHard = "yes", dzData = ATEPS_DZ, optimizer = "CSOLNP")
set.seed(1234)
ACEc <- umxConfint(ACEc, run = TRUE, level = .95, parm = c("top.A_std[1,1]", "top.C_std[1,1]", "top.E_std[1,1]"),
                  wipeExistingRequests = F,  optimizer = "SLSQP")
summary(ACEc, verbose = TRUE)

AE <- umxACE("AE", selDVs = "trust_sent_scaled", selCovs = c("male", "age"), sep = "_T", mzData = ATEPS_MZ, dzData = ATEPS_DZ, optimizer = "CSOLNP")
set.seed(1234)
AE = umxModify(AE, update = c("c_r1c1"), comparison = TRUE)
AE <- umxConfint(AE, run = TRUE, level = .95, parm = c("top.A_std[1,1]", "top.E_std[1,1]"), 
                  wipeExistingRequests = F, optimizer = "SLSQP")
summary(AE, verbose = TRUE)

CE <- umxACE("CE", selDVs = "trust_sent_scaled", selCovs = c("male", "age"), sep = "_T", mzData = ATEPS_MZ, dzData = ATEPS_DZ, optimizer = "CSOLNP")
set.seed(1234)
CE = umxModify(CE, update = c("a_r1c1"), comparison = TRUE)
CE <- umxConfint(CE, run = TRUE, level = .95, parm = c("top.C_std[1,1]", "top.E_std[1,1]"), 
                  wipeExistingRequests = F, optimizer = "SLSQP")
summary(CE, verbose = TRUE)

```



